home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / HENSA / MATHS / PLPLOT / PLPLOT.ZIP / examples / f77 / x09f.f < prev    next >
Encoding:
Text File  |  1994-05-27  |  6.4 KB  |  252 lines
  1. ! $Id: x09f.f,v 1.2 1994/05/26 19:34:23 mjl Exp $
  2. ! $Log: x09f.f,v $
  3. ! Revision 1.2  1994/05/26  19:34:23  mjl
  4. ! Inserted missing CVS Id and Log fields for all Fortran demos.  Comment
  5. ! character changed to "!" everywhere, to work well with font-lock in Lucid
  6. ! emacs (requires a small change to fortran-mode.el).
  7. !
  8. !
  9.       program example09
  10. !     =================
  11. !
  12. ! Demonstration of contour plotting      
  13.  
  14.       integer i,j,nptsx,nptsy
  15.       parameter (nptsx=41,nptsy=35)
  16.  
  17.       real z(nptsx,nptsy), w(nptsx,nptsy), clevel(11)
  18.       real xmin, xmax, ymin, ymax, x, y
  19.  
  20.       common /plplot/ tr(6)
  21.  
  22.       data clevel /-1.,-.8,-.6,-.4,-.2,0,.2,.4,.6,.8,1./
  23.       tr(1) = 0.05
  24.       tr(2) = 0.0
  25.       tr(3) = -1.0
  26.       tr(4) = 0.0
  27.       tr(5) = 0.05
  28.       tr(6) = -1.0
  29.  
  30.       x = 1.0
  31.       y = 1.0
  32.       xmin = tr(1) * (x-1) + tr(2) * (y-1) + tr(3)
  33.       ymin = tr(4) * (x-1) + tr(5) * (y-1) + tr(6)
  34.       x = nptsx
  35.       y = nptsy
  36.       xmax = tr(1) * (x-1) + tr(2) * (y-1) + tr(3)
  37.       ymax = tr(4) * (x-1) + tr(5) * (y-1) + tr(6)
  38.  
  39.       do 2 i=1,nptsx
  40.         xx = (i-1-(nptsx/2))/real(nptsx/2)
  41.         do 3 j=1,nptsy
  42.           yy = (j-1-(nptsy/2))/real(nptsy/2) - 1.0
  43.           z(i,j) = xx*xx - yy*yy
  44.           w(i,j) = 2*xx*yy
  45.     3   continue
  46.     2 continue
  47.  
  48.       call plinit()
  49.       call plenv(xmin,xmax,ymin,ymax,0,0)
  50.       call plcol(2)
  51.       call plcont(z,nptsx,nptsy,1,nptsx,1,nptsy,clevel,11)
  52.       call plstyl(1,1500,1500)
  53.       call plcol(3)
  54.       call plcont(w,nptsx,nptsy,1,nptsx,1,nptsy,clevel,11)
  55.       call plcol(1)
  56.       call pllab('X Coordinate', 'Y Coordinate', 
  57.      *           'Contour Plots of Saddle Points')
  58.  
  59.       call polar()
  60.  
  61.       call plend
  62.       end
  63.  
  64. !----------------------------------------------------------------------------!
  65. ! Routine for demonstrating use of transformation arrays in contour plots.
  66. ! Sorry for the formatting, as this has been through a preprocessor.
  67.  
  68.       subroutine polar()
  69.  
  70.       parameter (NCX=40, NCY=64, NPLT=200)
  71.       parameter (TWOPI=6.2831853071795864768)
  72.  
  73.       real z(NCX, NCY), ztmp(NCX, NCY+1)
  74.       real xg(NCX, NCY+1), yg(NCX, NCY+1), xtm(NPLT), ytm(NPLT)
  75.  
  76.       real clevel(40)
  77.       character*8 xopt, yopt
  78.  
  79.       nx = NCX
  80.       ny = NCY
  81.  
  82.       kx = 1
  83.       lx = nx
  84.       ky = 1
  85.       ly = ny
  86.  
  87. ! Set up r-theta grids
  88. ! Tack on extra cell in theta to handle periodicity.
  89.  
  90.       do 23000 i = 1, nx
  91.          r = i - 0.5
  92.          do 23002 j = 1, ny
  93.             theta = TWOPI/float(ny) * (j-0.5)
  94.             xg(i,j) = r * cos(theta)
  95.             yg(i,j) = r * sin(theta)
  96. 23002    continue
  97.          xg(i, ny+1) = xg(i, 1)
  98.          yg(i, ny+1) = yg(i, 1)
  99. 23000 continue
  100.       call a2mnmx(xg, nx, ny, xmin, xmax)
  101.       call a2mnmx(yg, nx, ny, ymin, ymax)
  102.       rmax = r
  103.       x0 = (xmin + xmax)/2.
  104.       y0 = (ymin + ymax)/2.
  105.  
  106. ! Potential inside a conducting cylinder (or sphere) by method of images.
  107. ! Charge 1 is placed at (d1, d1), with image charge at (d2, d2).
  108. ! Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2).
  109. ! Also put in smoothing term at small distances.
  110.  
  111.       eps = 2.
  112.  
  113.       q1 = 1.
  114.       d1 = r/4.
  115.  
  116.       q1i = - q1*r/d1
  117.       d1i = r**2/d1
  118.  
  119.       q2 = -1.
  120.       d2 = r/4.
  121.  
  122.       q2i = - q2*r/d2
  123.       d2i = r**2/d2
  124.  
  125.       do 23004 i = 1, nx
  126.          do 23006 j = 1, ny
  127.             div1 = sqrt((xg(i,j)-d1)**2 + (yg(i,j)-d1)**2 + eps**2)
  128.             div1i = sqrt((xg(i,j)-d1i)**2 + (yg(i,j)-d1i)**2 + eps**2)
  129.  
  130.             div2 = sqrt((xg(i,j)-d2)**2 + (yg(i,j)+d2)**2 + eps**2)
  131.             div2i = sqrt((xg(i,j)-d2i)**2 + (yg(i,j)+d2i)**2 + eps**2)
  132.  
  133.             z(i,j) = q1/div1 + q1i/div1i + q2/div2 + q2i/div2i
  134. 23006    continue
  135. 23004 continue
  136.  
  137. ! Tack on extra cell in theta to handle periodicity.
  138.  
  139.       do 23008 i = 1, nx
  140.          do 23010 j = 1, ny
  141.             ztmp(i,j) = z(i,j)
  142. 23010    continue
  143.          ztmp(i, ny+1) = z(i, 1)
  144. 23008 continue
  145.       call a2mnmx(z, nx, ny, zmin, zmax)
  146.  
  147. ! Set up contour levels.
  148.  
  149.       nlevel = 20
  150.       dz = abs(zmax - zmin)/float(nlevel)
  151.       do 23012 i = 1, nlevel
  152.          clevel(i) = zmin + (i-0.5)*dz
  153. 23012 continue
  154.  
  155. ! Split contours into two parts, z > 0, and z < 0.
  156. ! Dashed contours will be at levels 'ilevlt' through 'ilevlt+nlevlt'.
  157. ! Solid  contours will be at levels 'ilevgt' through 'ilevgt+nlevgt'.
  158.  
  159.       ilevlt = 1
  160.       nlevlt = 0
  161.       do 23014 i = 1, nlevel
  162.          if(clevel(i) .gt. 0.) go to 23015
  163.          nlevlt = nlevlt + 1
  164. 23014 continue
  165. 23015 continue
  166.       ilevgt = ilevlt + nlevlt
  167.       nlevgt = nlevel - nlevlt
  168.  
  169. ! Advance graphics frame and get ready to plot.
  170.  
  171.       ncollin = 11
  172.       ncolbox = 1
  173.       ncollab = 2
  174.  
  175.       call pladv(0)
  176.       call plcol(ncolbox)
  177.  
  178. ! Scale window to user coordinates.
  179. ! Make a bit larger so the boundary doesn't get clipped.
  180.  
  181.       eps = 0.05
  182.       xpmin = xmin - abs(xmin)*eps
  183.       xpmax = xmax + abs(xmax)*eps
  184.       ypmin = ymin - abs(ymin)*eps
  185.       ypmax = ymax + abs(ymax)*eps
  186.  
  187.       call plvpas(0.1, 0.9, 0.1, 0.9, 1.0)
  188.       call plwind(xpmin, xpmax, ypmin, ypmax)
  189.  
  190.       xopt = ' '
  191.       yopt = ' '
  192.       xtick = 0.
  193.       nxsub = 0
  194.       ytick = 0.
  195.       nysub = 0
  196.  
  197.       call plbox(xopt, xtick, nxsub, yopt, ytick, nysub)
  198.  
  199. ! Call plotter once for z < 0 (dashed), once for z > 0 (solid lines).
  200.  
  201.       call plcol(ncollin)
  202.       if(nlevlt .gt. 0) then
  203.          call pllsty(2)
  204.          call plcon2(ztmp, nx, ny+1, kx, lx, ky, ly+1, clevel(ilevlt), 
  205.      &nlevlt, xg, yg)
  206.       endif
  207.       if(nlevgt .gt. 0) then
  208.          call pllsty(1)
  209.          call plcon2(ztmp, nx, ny+1, kx, lx, ky, ly+1, clevel(ilevgt), 
  210.      &nlevgt, xg, yg)
  211.       endif
  212.  
  213. ! Draw boundary.
  214.  
  215.       do 23016 i = 1, NPLT
  216.          theta = (TWOPI)/(NPLT-1) * float(i-1)
  217.          xtm(i) = x0 + rmax * cos(theta)
  218.          ytm(i) = y0 + rmax * sin(theta)
  219. 23016 continue
  220.       call plcol(ncolbox)
  221.       call plline(NPLT, xtm, ytm)
  222.  
  223.       call plcol(ncollab)
  224.       call pllab(' ', ' ', 
  225.      &'Shielded potential of charges in a conducting sphere')
  226.  
  227.       return
  228.       end
  229.  
  230. !----------------------------------------------------------------------------!
  231. ! Subroutine a2mnmx
  232. !----------------------------------------------------------------------------!
  233. ! Minimum and the maximum elements of a 2-d array.
  234. !----------------------------------------------------------------------------!
  235.  
  236.       subroutine a2mnmx(f, nx, ny, fmin, fmax)
  237.  
  238.       integer nx, ny
  239.       real f(nx, ny), fmin, fmax
  240.  
  241.       fmax=f(1, 1)
  242.       fmin=fmax
  243.       do 23000 j=1, ny
  244.          do 23002 i=1, nx
  245.             fmax=max(fmax, f(i, j))
  246.             fmin=min(fmin, f(i, j))
  247. 23002    continue
  248. 23000 continue
  249.  
  250.       return 
  251.       end
  252.